from collections import defaultdict

def FastAlign(seq1, seq2):
    gapScore = -2
    ktup = 3
    cutoff = 2
    
    i = len(seq1)   #col - okresla kolumny
    j = len(seq2)   #row - okresla wiersze
    
    #slownik zawierajacy wszystkie gorace miejscela
    #klucze to numery przekatnych, wartosci to lista ze wspolrzednymi goracych miejsc
    hotDiagonals = defaultdict(list)
    
    #slownik currentScore'ow kazdej przekatnej (przekatna: currentScore)
    diagonalsScores = defaultdict(list)
    
    #tworzenie tablicy o j-wierszach oraz i-kolumnach
    matrix = [0]*j
    for row in xrange(len(matrix)):
        matrix[row] = [0]*i
        
    #pierwszy etap: przelatuje cala tablice w poszukiwaniu goracych miejsc
    for row in xrange(j):
        for col in xrange(i):
            #jesli sa takie same na obecnych pozycjach w obu sekwencjach 
            #to szukam goracych miejsc o dlugosci ktup
            if seq1[col] == seq2[row]:
                #funkcja sprawdza czy w tym miejscu jest gorace miejsce i jesli tak, to zaznacza je jedynkami
                markHotPlace(matrix, row, col, hotDiagonals, ktup, seq1, seq2)
                
    print "Dla sekwencji %s i %s \nznaleziono nastepujace gorace miejsca:" %(seq1, seq2)
    printMatrix(matrix)
    
    print "Gorace przekatne:"
    
    #kolejny etap algorytmu: szukam gapow i je zaznaczam
    for diagonalNum in hotDiagonals.keys():
        print "przekatna nr %i posiada gorace miejsca: %s" %(diagonalNum, hotDiagonals[diagonalNum]) 
        
        #oznaczanie przerw ma sens tylko gdy jest wiecej niz 1 gorace miejsce na przekatnej:
        if len(hotDiagonals[diagonalNum]) > 1:
            markGaps(diagonalNum, hotDiagonals[diagonalNum], matrix, ktup, gapScore)

    #znalazlem i oznaczylem gapy
    printMatrix(matrix)
    print "Wszystkieg gorace przekatne: ",hotDiagonals
    
    #obliczam scory dla przekatnych
    for diagonal in hotDiagonals.keys():
        diagonalScore = 0
        thisHotPlace = list(hotDiagonals[diagonal][0])
        lastHotPlace = list(hotDiagonals[diagonal][len(hotDiagonals[diagonal])-1])
            
        while thisHotPlace != lastHotPlace:
            diagonalScore += matrix[thisHotPlace[0]][thisHotPlace[1]]
            thisHotPlace[0] += 1
            thisHotPlace[1] += 1
        
        #w ^^^^^ tej petli pomijany jest ostatni goracy element, wiec teraz go dodamy    
        diagonalScore += matrix[lastHotPlace[0]][lastHotPlace[1]]
        diagonalsScores[diagonal] = diagonalScore
        
    print "Scory wszystkich przykatnych: ", diagonalsScores

    #sortujemy od najlepszej; (przekatna, ccore)
    diagonalsScores = sorted(diagonalsScores.iteritems(), key=lambda (k,v): (v,k), reverse=True)
    print "Scory najlepszej dziesiatki: ", diagonalsScores[:10]
    
    #rescoring 10 najlepszych przekatnych tym razem z tablica PAM lub BLOSUM
    for diagonal in diagonalsScores[:10]:
        for position in hotDiagonals[diagonal[0]]:
            matrix[position[0]][position[1]] = score(seq1[position[1]], seq2[position[0]])
    
    printMatrix(matrix)

    #uzupelniam tablice inits kolejnymi poduliniowieniami, a najlepsze to init1
    inits = []
    for diagonal in diagonalsScores[:10]:
        initScore = 0
        initPosition = []
        
        currentScore = 0
        currentPosition = [0,0]
        currentLenght = 0
        
        for step in hotDiagonals[diagonal[0]]:
            #sprawdzanie czy dalej jestem na ciaglym goracym miejscu w ramach jednej przekatnej
            if (step[0]-currentLenght) == currentPosition[0] and (step[1]-currentLenght) == currentPosition[1]:
                currentScore += matrix[step[0]][step[1]]
                
            else:
                #okazuje sie ze na tej przekatnej sa gapy, wiec tu je obslugujemy
                if currentScore > cutoff:
                    initScore = currentScore
                    initPosition = currentPosition
                
                currentScore = matrix[step[0]][step[1]]
                currentLenght = 0
                currentPosition = step
            
            #waruneczek taki ze znalezlismy cos lepszego od dotychczasowego init1
            if currentScore > cutoff:
                initScore = currentScore
                initPosition = currentPosition
            
            currentLenght += 1
          

        if initScore > cutoff:
            inits.append(((initPosition[0], initPosition[1]), currentLenght, initScore))
    
    #sortowanie tablicy initow wzgledem scora, aby init1 byl na poczatku:
    #TODO
    inits.sort(cmp)
    print "inity posortowane", inits
    
    #krok trzeci, buduje graf poduliniowien z dlugosciami sciezek pomiedzy nimi
    #klucz to wierzcholek1, wartosci to pary: wierzcholek2 i odleglosc(wierzcholek1,wierzcholek2)
    #graph = defaultdict(list)
    
    graph = {}

    for row in xrange(len(inits)):
        graphRoute = {}
        for col in xrange(len(inits)):
            vertex1 = inits[row]
            vertex2 = inits[col]
            print "vertexy", vertex1, vertex2
            
            #sprawdzam czy jest sciezka odpowiadajaca warunkom algorytmu
            if (vertex1[0][0]+vertex1[1]-1)<(vertex2[0][0]) and (vertex1[0][1]+vertex1[1]-1)<(vertex2[0][1]):
                print "wskoczylo"
                #jesli tak, to przypisuje jej wage zgodna z iloscia gapow
                graphRoute[inits[col]] = -gapScore*(vertex2[0][0]-vertex1[0][0]+vertex2[0][1]-vertex1[0][1])
                print graphRoute
        
        graph[inits[row]] = graphRoute        
   

    #bellman ford, do tego czasu jest wszystko ok
    allPaths = []
    allValues = []
    print "GRAF: ",graph
    
    for startVertex in graph.keys():
        d,p = bellman_ford(graph, startVertex)
        print "start ", startVertex
        print "\tp: ", p
        print "\td: ", d

       
        maximumDistance = float('-inf')
        minimumVertex = None
        path = []
        value = 0
        
        #szukanie minimum
        for v in d:
            if v != startVertex and d[v] != float('inf') and d[v] >= maximumDistance:
                maximumDistance = d[v]
                minimumVertex = v
                
            
        currentVertex = minimumVertex
        
        print "CURRENT VERTEX: ", currentVertex 
        
        if currentVertex == None:
            continue
        
        while startVertex != currentVertex:
            path.append(currentVertex)
            #wartosc poduliniowienia
            value += currentVertex[2]
            #wartosc gapa
            value -= d[currentVertex]
            currentVertex = p[currentVertex]
            
        path.append(currentVertex)
        path.reverse()
        allPaths.append(path)
        #zmiana znaku przy value:
        allValues.append(-value)
    
    #initn - best of the above ;)
    initnPath = None
    initnScore = 0
    
    for i in xrange(len(allPaths)):
        if allValues[i] > initnScore:
            initnScore = allValues[i]
            initnPath = allPaths[i]
       
    print "initN: ", initnPath, initnScore
    print "init1: ", inits[0]  
    
    #bierzemy te sekwencje:
    al1 = ""
    al2 = ""
    
    print initnPath[0]
    
    column = initnPath[0][0][1]
    row = initnPath[0][0][0]
    
    #odtworzenie calych poduliniowien
    for step in xrange(len(initnPath)):

        for i in xrange(initnPath[step][1]):
            al1 += seq1[column]
            al2 += seq2[row]
            
            row += 1
            column += 1
        
        if step < len(initnPath)-1:    
            next = (initnPath[step+1][0])
        else: 
            break
        
        while True:
            if (row, column) == next:
                break
            elif row < next[0]:
                row += 1
                al1 += "-"
                al2 += seq2[row]
                
            elif column < next[1]:
                column += 1
                al1 += seq1[column]
                al2 += "-"
          
    print "KONCOWE LOKALNE DOPASOWANIE SEKWENCJI INITN:"      
    print al1
    print al2
    
    #smith&waterman
    
    #DynAlign(al1, al2, -2) 
        
    #i po belmanie
    
def initialize(graph, source):
    d = {}
    p = {}
    for node in graph:
        d[node] = float('Inf')
        p[node] = None
    d[source] = 0
    
    return d, p

def relax(u, v, graph, d, p):
    if d[v] > d[u] + graph[u][v]:
        d[v]  = d[u] + graph[u][v]
        p[v] = u

def bellman_ford(graph, source):
    d, p = initialize(graph, source)
    for i in range(len(graph)-1):
        for u in graph:
            for v in graph[u]:
                relax(u, v, graph, d, p)
    for u in graph:
        for v in graph[u]:
            assert d[v] <= d[u] + graph[u][v]
    return d, p
 
    #koniec belmana, dalej jest juz wszystko ok

#TODO zmienic kiedys oceny na tablice substytucji
def score(x, y):
    if x == y:
        return 4
    else:
        return -1

#funkcja sortujaca liste krotek wzgledem drugiego elementu (scora), malejaco 
def cmp(a,b):
    if a[2] > b[2]:
        return -1
    elif a[2] < b[2]:
        return 1
    else:
        return 0

#ustala kary za przerwy pomiedzy goracymi miejscami na tej samej przekatnej
def markGaps(diagonalNum, hotPlaces, matrix, ktup, gapScore):
    # gapScore = -2
    
    for placeNumber in xrange(len(hotPlaces)):
        thisRow = hotPlaces[placeNumber][0]
        thisCol = hotPlaces[placeNumber][1]
        
        #sprawdzam czy nie wyszedlem poza tablice
        if placeNumber+1 == len(hotPlaces) or thisRow+1 == len(matrix) or thisCol+1 == len(matrix[0]):
            return
        
        while (thisRow+1, thisCol+1) != hotPlaces[placeNumber+1]:
            print "znalazlem gap w: ", (thisRow+1, thisCol+1)
            
            matrix[thisRow+1][thisCol+1] = gapScore
            thisRow += 1
            thisCol += 1
    
#sorawdza czy w matrix[row][col] jest gorace miejsce o dlugosci ktup
#jesli tak, to gorace miejsce wystepuje w tablicy jako 1
def markHotPlace(matrix, row, col, hotDiagonals, ktup, seq1, seq2):
    #score dla goracego miejsca
    hotPlaceScore = 1
    
    for step in xrange(ktup):
        #sprawdzam czy nie wyszlismy poza tablice
        if (row+step) == len(matrix) or (col+step) == len(matrix[0]):
            return
        #nie ma takich samych literek w seq1 i seq2 => brak goracych miejsc
        elif seq1[col+step] != seq2[row+step]:
            return
    
    #jesli przeszedlem przez poprzednie etapy caly
    #to znaczy ze znalazlem gorace miejsca, teraz przypisuje im score
    for step in xrange(ktup):
        matrix[row+step][col+step] = hotPlaceScore
        #zapisuje gorace miejsce
        if (row+step, col+step) not in hotDiagonals[col-row]:
            hotDiagonals[col-row].append((row+step, col+step))
    
    
        
#drukowanko
def printMatrix(matrix):
    for row in xrange(len(matrix)):
        print matrix[row]
        
        
#bez komentarza...
def DynAlign(seq1, seq2, gap):
    # tu chyba tez wszystko jasne
    m = len(seq1)+1
    n = len(seq2)+1
    
    # tu tworze tablice dwuwymiarowa m*n
    # musi byc o jeden wieksza od rozmiarow sekwencji aby w pierwszej kolumnie i wierszu byly zera
    # zapewniajace zakonczenie odtwarzania dopasowania gdy sie na nie trafi
    scoreTable = [0]*m  
    for i in range(len(scoreTable)):
        scoreTable[i] = [0]*n
    
    maxValue = 0
    maxValuePosition = (0,0)
    
    # glowna petla programu wypelniajaca tablice podobienstw zgodnie z funcja score i karami za przerwe
    for row in range(1, m):
        for col in range(1, n):
            # poprzednia wartosc po skosie + punkty za obecne elementy 
            v1 = scoreTable[row-1][col-1] + score(seq1[row-1], seq2[col-1])
            # poprzednia wartosc z gory + kara za przerwe
            v2 = scoreTable[row][col-1] + gap
            # poprzednia wartosc z lewej + kara za przerwe
            v3 = scoreTable[row-1][col] + gap
            # tego chyba nie trzeba tlumaczyc...
            scoreTable[row][col] = max(v1, v2, v3, 0)
            #wybieranie maksymalnej wartosci z calej tablicy
            if scoreTable[row][col] >= maxValue:
                maxValue = scoreTable[row][col]
                maxValuePosition = (row, col)             
    
    row = maxValuePosition[0]
    col = maxValuePosition[1]
    align1 = ""
    align2 = ""
    
    # petla odtwarzajaca dopasowane ciagi
    # w danym punkcie tablicy sprawdzam czy przyszedlem do niego odpowiednio:
    # czy z gory + gap
    # czy z lewej + gap
    # else: po skosie + punkty
    while scoreTable[row][col] > 0:
        #sprawdzam czy przyszedlem z gory
        if scoreTable[row][col] == scoreTable[row-1][col] + gap:
            align1 = seq1[row-1] + align1
            align2 = "-" + align2
            row -= 1
        
        #tu: czy przyszedlem z lewej
        elif scoreTable[row][col] == scoreTable[row][col-1] + gap:
            align1 = "-" + align1
            align2 = seq2[col-1] + align2
            col -= 1
        
        #w ostatecznosci przyszedlem po skosie    
        else:
            align1 = seq1[row-1] + align1
            align2 = seq2[col-1] + align2
            row -= 1
            col -= 1
    




seq1 = "TTTGTTTTGCTGTGTCTGATTGATTGATTCATTTTTAAGACACTGACTATCCCATA"
seq2 = "CTCTGTTCAGCCCCACGGTTAGGCTAA"

#FastAlign(seq1, seq2)
#FastAlign("TTTGTTTTGCTGTGTCTGATTGATTGATTCATTTTTAAGACACTGACTATCCCATA", "CTCTGTTCAGCCCCACGGTTAGGCTAAGATCTGTTCTACTTTCTCCCTAGGGCTTAATTGCCAGACTGACTTACCGGTGCCTCCTAGCCTTTGTTCTGCCTCTCACGCTGACATCAGGGAAGCAGTGATGTAGAATAAAGCAGGGGTG")
